MODULE ModCoolOI

!!$ This module contains the subroutine for calculating the coolig due to OI.
!!$ Programming conventions are from Nyhoff and Leestma, _Fortran 90 for Scientists and Engineers_, 1997
!!$
!!$ STH: Samuel Harrold

!!$ Declare modules to USE
  USE ModParms
  USE ModFuncs
  USE ModAbund

!!$ Do not assume undeclared variables
  IMPLICIT NONE



CONTAINS



!!$ CalcCoolOI
!!$ Calculate the cooling due to the magnetic dipole and electric quadrapole fine-structure forbidden transition [OI] 2P0 -> 2P2 (E2), 2P0 -> 2P1 (M1), 2P1 -> 2P2 (M1). Energy levels from highest to lowest: 2P0, 2P1, 2P2. The cooling function follows the treatment described in Tielens and Hollenbach 1985 (TH85), the collisonal coefficients also from TH85, and the solution for level populations from White 1961. Transitions rates denoted as rate from upper to lower = Ruprlwr (e.g. R02)
!!$ Accepts: particle number density, temperature, mean intensity of radiation field at each transition with respect to cm^-1, column density in line of sight to the star
!!$ Returns: cooling of the ga in erg s^-1 (gas particle)^-1, abundance of OI as TEST
!!$ TODO: calculate cooling by arrays rather than by element

  SUBROUTINE CalcCoolOI(ngas, Tempgas, Jw02, Jw01, Jw12, Ntostar, CoolOI)

!!$ Declare variables
!!$

!!$ Specify input and output
!!$

!!$ ngas       = density of gas in particles cm^-3
!!$ Tempgas    = temperature of gas in K
!!$ Jw02/01/12 = mean intensity of radiation field in erg/(s cm^2 sr cm^-1)
!!$ for each transition. w for omega as wavenumber
!!$ Ntostar    = column density of gas in line of sight to star in particles cm^-2
!!$ CoolOI     = cooling due to emission from OI in units acceptable by pisco 
!!$ in erg s^-1 (gas particle)^-1
    REAL, INTENT(IN)  :: ngas , Tempgas, Jw02, Jw01, Jw12, Ntostar
    REAL, INTENT(OUT) :: CoolOI

!!$ Specify internal variables
!!$

!!$ Data for [OI] at 44, 63, 145 microns from NIST Atomic Spectra Database
!!$ A02/01/12  = radiative de-excitation coefficients for each transition in s^-1
!!$ g0/1/2     = statistical weight g for each level (g = 2J + 1)
!!$ nu02/01/12 = frequency of each OI transition from microns into Hz
    REAL :: &
         A02 = 1.34E-10, &
         A01 = 1.75E-5, &
         A12 = 8.91E-5, &
         g0  = 1., &
         g1  = 3., &
         g2  = 5., &
         nu02 = c / (44.0573  * 1.E-4), &
         nu01 = c / (145.5350 * 1.E-4), &
         nu12 = c / (63.1852  * 1.E-4)
    
!!$ Variables for physical conditions
!!$ Notation follows that of Tielens and Hollenbach 1985 (TH85).
!!$ Jnu02/01/12 = mean intensity of radiation field wrt Hz, units erg/(s cm^2 sr Hz), for each transition
!!$ Abund(OI/p/H/H2) = abundances of OI/p/H/H2 relative to all gas particles determined by column density to the star
!!$ n(OI/p/H/H2)     = number density of OI (all levels)/p/H/H2
!!$ Besc             = escape probability
!!$ Pnu02/10/12      = background radiation field wrt Hz for each transition
!!$ Qnu02/01/12      = population mode of background radiation field wrt Hz for each transition
!!$ gamma(ij in combo{0,1,2})(p/H/H2) = collisional excitation/de-excitation coefficients for all colliding species. p from Tielens 2005. H from Launay, Roeff 1977b via TH85. H2 from Launay and Roeff 1977b via TH85.
!!$ C(ij in combo{0,1,2}) = composite collisional excitation/de-excitation coefficient for al transitions (atom s^-1 cm^-3)
!!$ R(ij in combo{0,1,2}) = rate for all transition (atom s^-1 cm^-3)
!!$ CoR(00,11,22) = cofactor for transtion rate matrix from equations of statistical equilibrium to solve for population densities (see White 1961)
!!$ lambda        = constant linking cofactor from conservation of particles (White 1961)
!!$ n(0/1/2)      = population density of state 0/1/2
!!$ Snu(02/01/12) = source function wrt Hz for all transitions
!!$ L(02/01/12) = cooling from TH85, Eqn B1 in units of erg s^-1 cm^-3 for each transition
!!$ TODO: include cofactors
    REAL :: &
         Jnu02, Jnu01, Jnu12, &
         AbundOI, Abundp, AbundH, AbundH2, &
         nOI, np, nH, nH2, &
         Besc, &
         Pnu02, Pnu01, Pnu12, &
         Qnu02, Qnu01, Qnu12, &
         gamma02p,  gamma01p,  gamma12p, &
         gamma20p,  gamma10p,  gamma21p, &
         gamma02H,  gamma01H,  gamma12H, &
         gamma20H,  gamma10H,  gamma21H, &
         gamma02H2, gamma01H2, gamma12H2, &
         gamma20H2, gamma10H2, gamma21H2, &
         C02, C01, C12, &
         C20, C10, C21, &
         R02, R01, R12, &
         R20, R10, R21, &
         CoR00, CoR11, CoR22, &
         lambda, &
         n0, n1, n2, &
         Snu02, Snu01, Snu12, &
         L02, L01, L12

!!$ Convert mean intensity of radiation field
!!$ Jnu = mean intensity of radiation field converted to erg/(s cm^2 sr Hz)
!!$ using nu = c*w => dnu = c*dw and Jnu is defined as a derivative
!!$ wrt nu, Jnu = dJ/dnu = dJ/dw * dw/dnu = Jw/c
    Jnu02 = Jw02 / c
    Jnu01 = Jw01 / c
    Jnu12 = Jw12 / c

!!$ Estimate number densities of species
!!$ number density of species = abundance * number density of gas
!!$ Note: assuming that majority of positive ions collide with H and transfer charge so that p abundance = e abundance (from Lacy). p more important colliders with atoms from Tielens 2005, Table 2.7
    CALL CalcAbund(Ntostar = Ntostar, &
         AbundOI = AbundOI, Abunde = Abundp, &
         AbundH = AbundH, AbundH2 = AbundH2)
    nOI = ngas * AbundOI
    np  = ngas * Abundp
    nH  = ngas * AbundH
    nH2 = ngas * AbundH2

!!$ Calculate escape probability, Besc
!!$ see de Jong, Dalgarno, and Boland 1980 for plane parallel approximation
!!$ Assume escape probability is 0.5 from face of disk
!!$ TODO: make tauIR an input
    Besc = 0.5

!!$ Next calculations follow Tielens and Hollenbach, 1985
!!$

!!$ Calculate background radiation field for each transition wrt Hz in erg / (s cm^2 sr Hz)
    Pnu02 = Jnu02
    Pnu01 = Jnu01
    Pnu12 = Jnu12

!!$ Calculate population mode of radiation field for each transition (unitless)
    Qnu02 = (c**2 / (2*h*nu02**3)) * Pnu02
    Qnu01 = (c**2 / (2*h*nu01**3)) * Pnu01
    Qnu12 = (c**2 / (2*h*nu12**3)) * Pnu12

!!$ Calculate collisional coefficient for excitation and deexcitation for each transition and species. Listed from most to least important
!!$ Collision coefficient gamma p from Tielens 2005
    gamma02p = 1.4E-8
    gamma20p = (g0/g2) * gamma02p * EXP(-h*nu02 / (k*Tempgas))
    gamma01p = 5.0E-9
    gamma10p = (g0/g1) * gamma01p * EXP(-h*nu01 / (k*Tempgas))
    gamma12p = 1.4E-8
    gamma21p = (g1/g2) * gamma12p * EXP(-h*nu12 / (k*Tempgas))

!!$ gamma H from Launay and Roeff 1977b via Tielens, Hollenbach 1985
    gamma02H = 1.1E-12 * Tempgas**0.8
    gamma20H = (g0/g2) * gamma02H * EXP(-h*nu02 / (k*Tempgas))
    gamma01H = 1.5E-10 * Tempgas**0.44
    gamma10H = (g0/g1) * gamma01H * EXP(-h*nu01 / (k*Tempgas))
    gamma12H = 4.2E-12 * Tempgas**0.67
    gamma21H = (g1/g2) * gamma12H * EXP(-h*nu12 / (k*Tempgas))

!!$ gamma H2 from Launay and Roeff 1977b via Tielens, Hollenbach 1985
    gamma02H2 = gamma02H
    gamma20H2 = gamma20H
    gamma01H2 = gamma01H
    gamma10H2 = gamma10H
    gamma12H2 = gamma12H
    gamma21H2 = gamma21H

!!$ Calculate collisional excitation/de-excitation rates. Terms are listed from most to least important.
    C02 = gamma02p*np + gamma02H*nH + gamma02H2*nH2
    C20 = gamma20p*np + gamma20H*nH + gamma20H2*nH2
    C01 = gamma01p*np + gamma01H*nH + gamma01H2*nH2
    C10 = gamma10p*np + gamma10H*nH + gamma10H2*nH2
    C12 = gamma12p*np + gamma12H*nH + gamma12H2*nH2
    C21 = gamma21p*np + gamma21H*nH + gamma21H2*nH2

!!$ Calculate excitation/de-excitation rates, units transition/(s cm^3)
    R02  = A02*Besc*(1. + Qnu02)  + C02
    R20  = (g0/g2)*A02*Besc*Qnu02 + C20    
    R01  = A01*Besc*(1. + Qnu01)  + C01
    R10  = (g0/g2)*A01*Besc*Qnu01 + C10
    R12  = A12*Besc*(1. + Qnu12)  + C12
    R21  = (g0/g2)*A12*Besc*Qnu12 + C21        

!!$ Calculate cofactors from rate equations for statistical equilibrium following Rosseland 1926, 1936 via White 1961.
!!$ Note from White 1964 notation, Pjm is cofactor of column j, row m.
    CoR00 = R10*R20 + R12*R20 + R10*R21
    CoR11 = R01*R20 + R01*R21 + R02*R21
    CoR22 = R02*R10 + R01*R12 + R02*R12

!!$ Calculate coefficient from conservation equation following Rosseland 1926, 1936 viaWhite 1961
    lambda = nOI / (CoR00 + CoR11 + CoR22)

!!$ Calculate population level densities, units atoms/cm^3
    n0 = lambda * CoR00
    n1 = lambda * CoR11
    n2 = lambda * CoR22

!!$ Calculate source function wrt Hz for each transition, units erg/(s cm^2 sr Hz)
    Snu02 = (2.*h*nu02**3/c**2) * (( (g0*n2)/(g2*n0) ) - 1.)**(-1)
    Snu01 = (2.*h*nu01**3/c**2) * (( (g0*n1)/(g1*n0) ) - 1.)**(-1)
    Snu12 = (2.*h*nu12**3/c**2) * (( (g1*n2)/(g2*n1) ) - 1.)**(-1)

!!$ Calculate cooling wrt Hz for each transition, units erg/(s cm^3)
    L02 = n0*A02*h*nu02*Besc*(1. - (Pnu02/Snu02))
    L01 = n0*A01*h*nu01*Besc*(1. - (Pnu01/Snu01))
    L12 = n1*A12*h*nu12*Besc*(1. - (Pnu12/Snu12))

!!$ Calcuate total cooling due to OI, units erg s^-1 (gas particle)^-1
    CoolOI = (L02 + L01 + L12) / ngas

  END SUBROUTINE CalcCoolOI



END MODULE ModCoolOI
